% rsts_start.M
%
% finds starting points for MCMC chains 
%
% Created: 12/11/2007 by Taeyoung Doh 

%--------------------------------------------------------------------
% housekeeping
%--------------------------------------------------------------------
%clear all;
close all;
diary off;


%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec2;

% check directories
retdir = chkdir(path_.root);
retdir = chkdir(path_.post);
retdir = chkdir(path_.result);


% parameter names
pnames = feval(fnames_.pnames,msel);

%--------------------------------------------------------------------
% prior specification
%--------------------------------------------------------------------
global prior_;
prior_ = feval(fnames_.prior,msel);

global mprior; 
%--------------------------------------------------------------------
% data and parameters
%--------------------------------------------------------------------

% load data
loaddata;

% paremeters
mhrun=0; 
path_.mhrun = [path_.para 'mhrun' num2str(mhrun) '/'];

%  Load Prior draws 
nblock=10; 

modsize=1000; 
llikstart=zeros(nblock,1); 
nm=100; 

% Looking for a starting point for posterior density maximization 
% Pick up the parameter draw with the highest likelihood in each block of
% prior draws 
paraast=[];
llikast=[];  
for iblock=1:nblock 
llikpara=[];
pfile=[path_.mhrun 'mhdraw' num2str(iblock,'%04.0f') '.mat']; 
paras=load(pfile); 
parast=paras.parasim; 
index =paras.retcode==0;
parast=parast(index,:); 
nparasim=length(parast); 
 for j=1:floor(nparasim/nm):nparasim
     
           llikpara(j,1) = tloglik4(parast(j,:)',data);    
      
 end  
 % print output
fmt='%15.8f ';		% printing format

disp(' ');
disp(sprintf('block: %2.0f / %2.0f        simulation steps: %5.0f / %5.0f',iblock,nblock,j,nparasim));
disp(sprintf('-----------------------------------------------------'));
disp(sprintf(['maximum likelihood in this block  : ' fmt],max(llikpara)));
[B,l]=sort(llikpara,1,'descend');    
paraast = [paraast ; parast(l(1),:)];
llikstart(iblock)=max(llikpara); 
end


% save output for later uses
sfile = [path_.post 'start.mat'];
save(sfile,'paraast','llikstart');
% global variables for parameter specification
global paraspec_;
paraspec_ = feval(fnames_.prtr,msel);


% initial parameters to maximize
[C,stind]=max(llikstart); 
para  = paraast(stind,:)'; 
xinit = invtrans(para);			% 'model' to 'max'
Hinit = eye(length(para))*1e-8;	% initial Hessian, no information
npara = length(para);

%--------------------------------------------------------------------
% optimization
%--------------------------------------------------------------------

% set optimization routine environment
nit  = 1000;		% maximum number of iteration 
crit = 1E-12;		% tolerance 

% We need to define 'cost' function to anneal
% and it is preferrably a function in 'max' parameters 
% to be used with CSMINWEL.
% See FNLINCOST.M for details.
% We use global variable DATA to make it simple 
% to deal with 'data' argument in CSMIN procedure.
% (You don't need to pass it as an additional argument by doing this.)

% Run CSMINWEL
%
% INPUT
%	'fnlincost' : cost function, FNLINCOST.M
%	xinit        : initial parameter values
%   Hinit        : initial Hessian
%   crit         : tolerance
%   nit          : maximum number of iteration
%
% OUTPUT
%   pmin         : minimized cost function value
%   xstar        : minimizer
%   pmgrad       : gradient
%   pmhess       : Hessian
%   itct         : iteration count
%   fcount       :
%   retcode      : return code, see CSMINIT.M for detail

[pmin,xstar,pmgrad,pmhess,itct,fcount,retcode] = csminwel('fncost',xinit,Hinit,[],crit,nit,data);


%--------------------------------------------------------------------
% posterior mode
%--------------------------------------------------------------------
% posterior mode in 'model' parameters
pmod = trans(xstar);

% posterior value at the mode
pmax = rsfnmh(pmod,data);
plik = tloglik4(pmod,data); 
ppri = logprior(pmod);

% Compute covariance matrix to be used for proposing step 
% in Metropolis-Hastings algorithm. It is in 'max' parameters.
[sigmult,penalt] = mhcov(pmod,data);

% compute std of mode estimates
pmod_std = sqrt(diag(sigmult*sigmult'));

%--------------------------------------------------------------------
% save results and print
%--------------------------------------------------------------------

% save output for later uses
sfile = [path_.post 'pm.mat'];
save(sfile,'pmod','pmax','plik','ppri','penalt','sigmult','pmod_std');


% print output
fmt='%15.8f ';		% printing format

% header
disp(' ');
disp('==================================================================')

% posterior mode
disp(' ');
disp('Parameter  Fixed           Mode           StDev           ');
for i=1:length(para)
	disp(sprintf(['%s      %d'  repmat(fmt,1,2) ],pnames(i,:),prior_.mask(i),pmod(i),pmod_std(i)));
end
disp(' ');
disp(sprintf(['Log Prior      : ' fmt],ppri));
disp(sprintf(['Log Likelihood : ' fmt],plik));
disp(sprintf(['Log Posterior  : ' fmt],pmax));


